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Abstract. 

Performing fully general relativistic simulations taking account of microphysical 
processes (e.g., weak interactions and neutrino cooling) is one of long standing problems 
in numerical relativity. One of main difficulties in implementation of weak interactions 
in the general relativistic framework lies on the fact that the characteristic timescale 
of weak interaction processes (the WP timescale, t wp ~ |Y^,/Y e |) in hot dense matters 
is much shorter than the dynamical timescale (idyn)- Numerically this means that 
stiff source terms appears in the equations so that an implicit scheme is in general 
necessary to stably solve the relevant equations. Otherwise a very short timestep 
(At < t wp <C idyn) will be required to solve them explicitly, which is unrealistic in 
the present computational resources. Furthermore, in the relativistic framework, the 
Lorentz factor is coupled with the rest mass density and the energy density. The 
specific enthalpy is also coupled with the momentum. Due to these couplings, it is very 
complicated to recover the primitive variables and the Lorentz factor from conserved 
quantities. Consequently, it is very difficult to solve the equations implicitly in the fully 
general relativistic framework. At the current status, no implicit procedure have been 
proposed except for the case of the spherical symmetry. Therefore, an approximate, 
explicit procedure is developed in the fully general relativistic framework in this paper 
as an first implementation of the microphysics toward a more realistic sophisticated 
model. The procedure is based on the so-called neutrino leakage schemes which is 
based on the property that the characteristic timescale in which neutrinos leak out 
of the system (the leakage timescale, £i ea k) is much longer than the WP timescale. 
In the previous leakage schemes, however, the problems of the stiff source terms are 
avoided in an artificial manner. In this paper, I present a detailed neutrino leakage 
scheme and a simple and stable method for solving the equations explicitly in the fully 
general relativistic framework. The drawback of the artificial treatment of the stiff 
source terms is improved. I also perform a test simulation to check the validity of the 
present method, showing that it works fairly well. 
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1. Introduction 

1.1. Motivation 

Numerical relativity is the unique and powerful tool to explore dynamical phenomena 
in which strong gravity plays important roles. Stellar core collapse and mergers of 
compact star binaries are among the most important and interesting events in the 
field. In theoretical view points, performing simulations of these phenomena is one 
of challenging problems because a rich diversity of physics has to be taken into account. 
All four known forces of nature are involved and play important roles during the collapse. 
General relativistic gravity plays essential roles in formation of a black hole. Note also 
that general relativity may play important roles in supernova explosion as previous 
pioneering works [U |2] showed. The weak interactions and emission of neutrinos govern 
energy and lepton-number losses, and hence driving the thermal and chemical evolutions 
of the system. The strong interactions determine ingredients and properties of dense 
matters. Strong magnetic fields, if they are present, may modify the dynamics. 

There is a long list of studies which explore these phenomena in the framework of 
numerical relativity (see [3] and references therein for resent simulations of stellar core 
collapse, and see jH [5] and references therein for those of compact binary merger). In 
most of the previous studies, however, treatments of microphysics are very simplified 
and more sophisticated studies are necessary. 

Furthermore, recent observations [6] [7J [8j [9j [TO], [TT] [12] have discovered the 
spectroscopic connections between several supernovae and long gamma-ray bursts 
(GRBs), clarifying that at least some of long GRBs are associated with the collapse 
of massive stars. Also, there are theoretical models that a short GRB occurs as a result 
of a binary neutron star merger [131 E3] • The relevant process of the energy deposition 
to form a GRB fireball may be pair annihilation of neutrinos emitted from a hot massive 
disk around a black hole formed after the collapse or the merger. These also enhance 
the importance of exploring stellar core collapse and coalescence of compact star binary 
in full general relativity taking account of microphysical processes. 

Gravitational wave astronomy will start in this decade. The first generation of 
ground-based interferometric detectors (LIGO [15J, VIRGO [16J, GEO600 [17]) are 
now in the scientific search for gravitational waves. To obtain physically valuable 
information from these observations, it is necessary to connect the observed data and 
the physics behind it. For this purpose, performing numerical simulation is the unique 
approach. However, accurate predictions of gravitational waveforms are still hampered 
by the facts that reliable estimates of waveforms require a general relativistic treatment 
[18] [19], and that appropriate treatments of microphysics such as a nuclear equation 
of state (EOS), the electron capture, and neutrino emissions and transfers. General 
relativistic simulations including microphysics are required to make accurate predictions 
of gravitational waveforms. 

As described above, to perform multidimensional simulations in the frame work of 
numerical relativity implementing microphysics is currently one of the most important 
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subjects in theoretical astrophysics. In spherical symmetry, fully general relativistic 
simulations [201 12U 122] of stellar core collapse have been performed in so-called state- 
of-the-art manners, namely, employing realistic equations of state, taking account of 
relevant microphysics, and solving the Boltzmann equation for the transfer of neutrinos 
(see also [23] ). In the multidimensional case, by contrast, there are few studies in the 
framework of general relativity. Recentlygeneral relativistic simulations implementing 
a realistic EOS and the electron capture were performed [2H [25] . In their calculation, 
however, the electron capture rate is not calculated in a self-consistent manner. Instead, 
they adopted a simplified prescription proposed in |26j, which is based on results of 
spherically symmetric simulations. It is not clear whether this treatment is justified for 
non-spherical collapse and the multidimensional phenomena. More importantly, they 
did not take account of neutrino cooling. 

Recently, I have made a fully general relativistic code with microphysics for the first 
time [2Z]. Since it is currently impossible to fully solve the multidimensional neutrino 
transfer equations in the framework of full general relativity because of restrictions of 
computational resources, it will be reasonable to adopt an approximated treatment of 
neutrino cooling. In that work, a general relativistic version of the so-called neutrino 
leakage schemes is developed. 

1.2. Neutrino leakage scheme 

The leakage schemes [28J |29j EH ED [32] as an approximate method for the neutrino 
cooling has a well-established history (e.g. [31J). The basic concept of the original 
neutrino leakage schemes [28, 29J is to treat the following two regions in the system 
separately: one region is where the diffusion timescale of neutrinos is longer than the 
dynamical timescale, and hence neutrinos are 'trapped' (neutrino-trapped region); the 
other region is where the diffusion timescale is shorter than the dynamical timescale, and 
hence neutrinos stream out freely out of the system (free-streaming region). The idea 
of treating the diffusion region separately has been applied to more advanced methods 
for the neutrino transfer (see e.g., |33j and references therein). 

Then, electron neutrinos and anti-neutrinos in the neutrino-trapped region are 
assumed to be in /9-equilibrium state. The net local rates of lepton-number and energy 
exchange with matters are set to be zero in the neutrino-trapped region. To treat 
diffusive 'leakage' of neutrinos out of the neutrino-trapped region, phenomenological 
source terms based on the diffusion theory are introduced [28l [29] . In the free-streaming 
region, on the other hand, it is assumed that neutrinos escape from the system without 
interacting with matter. Therefore, neutrinos carry the lepton number and the energy 
according to the local weak interaction rates. The neutrino fractions are not solved in 
the original version of the leakage scheme: Only the total lepton fraction is necessary in 
the free-streaming region and the neutrino fractions are set to zero in the free-streaming 
region. Note that there is a sharp discontinuity between the two regions. Consequently, 
thermodynamical quantities, in particular those of neutrinos and the electron fraction, 
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are also discontinuous at the boundary. 

The transfer of neutrinos are not solved in the leakage schemes. Therefore, they 
cannot treat non-local interactions among the neutrinos and matters; for example, the 
so-called neutrino heating [31] and the neutrino pair annihilation cannot be treated in the 
leakage scheme. Nevertheless, I consider a detailed general relativistic leakage scheme 
presented in this paper to be an important step towards more reliable and sophisticated 
models, since the simulated physical timescales in the case of compact binary mergers 
will be order of 10 ms and neutrino transfer is expected to be unimportant [35], and 
since the neutrino heating would be not very important in the case of prompt black hole 
formation. 

Usually, the boundary between the neutrino-trapped and free-streaming regions is 
given by hand as a single 'neutrino-trapping' density (ptrap) m the previous simulations 
of stellar core collapse [281 1291 132]. In fact, however, the location at which the neutrino 
trapping occurs depends strongly on the neutrino energies (e u ), and hence, there are 
different neutrino-trapping densities for different neutrino energies. The neutrino- 
trapping densities depend strongly on the neutrino energies as p trap oc e~ 3 [36]. This 
implies that neutrinos with lowest energy leave their corresponding neutrino-trapping 
region first, and neutrinos with higher energy are emitted later. 

In the previous leakage schemes [28l 1291 132] . on the other hand, all neutrinos are 
emitted in one moment irrespective of their energy. Consequently in the case of the 
so-called neutrino burst emission (e.g., [36]), for example, the duration in which the 
neutrinos are emitted is shortened and the peak luminosity at the burst is overestimated 
in the previous leakage schemes [29j [27]. The dependence of the neutrino-trapping 
densities and neutrino diffusion rates on the neutrino energies are approximately taken 
into account in the recent simulations of binary neutron star mergers J37J [35] . However 
the lepton- number conservation equations for neutrinos are not solved [37] . 

Recently, a numerical code based on a relativistic extension of the leakage schemes 
was developed in [27], where not the region of the system but the energy momentum 
tensor of neutrinos are decomposed into two parts; 'trapped-neutrino' and 'streaming- 
neutrino' parts. However the source terms of hydrodynamic and the lepton-number- 
conservation equations are determined using the single neutrino-trapping density as 
in the case of the previous leakage schemes. More recently, Liebendorfer et al. [38] 
proposed a scheme, which they call the isotropic diffusion source approximation, where 
the neutrino distribution function is decomposed into an isotropic distribution function 
of trapped neutrinos and a distribution function of streaming neutrinos. 

The present work is based on the previous studies described above. The framework 
of general relativistic extension of leakage scheme is based on my previous study in 
[27]. The treatment of neutrino diffusion rates is based on the recent work by Rosswog 
and Liebendorfer [35] where the neutrino-energy dependences are taken into account. 
Thus the remaining main problem to implement the relevant microphysics is that 
straightforward explicit scheme cannot be adopted to solve the equations [10] since 
the characteristic timescale of weak interaction processes (t wp ~ \Y e /Y e \, hereafter the 
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WP timescale) is much shorter than the dynamical timescale (idyn) m hot dense regions, 
as described in Sec. [2j Note that the WP timescale is different from the so-called weak 
timescale. In this paper I present a simple and stable method in which the equations are 
solved explicitly in the dynamical timescale in the fully general relativistic framework. 

The paper is organized as follows. First, main difficulties of implementation of weak 
interactions and neutrino cooling in full general relativity compared to implementation 
of them in Newtonian framework are briefly summarized in Sec. |2j Then, framework 
of the implementation of the microphysics is described in detail in Sec. El Some details 
of microphysics and numerics are described in Sec. HJ although GR leakage framework 
is independent of specific implementations of microphysics. In Sec. El results of a test 
simulation is briefly described illustrating good ability of the present implementation. 
Section |6] is devoted to the summary and discussions. Throughout the paper, the 
geometrical unit c = G = 1 is used otherwise stated. 

2. Difficulties of implementation of weak interactions and neutrino cooling 
in full general relativity 

Since the characteristic timescale of weak interaction precesses (the WP timescale 
tw P ~ | Ye/y^l) is much shorter than the dynamical timescale tdyn in hot dense matters 
[40. 35J, the numerical treatment of the weak interactions cannot be explicit, as noted 
in the previous pioneering work by Bruenn [10] ■ Otherwise a very short timestep (At < 
t wp <C ^dyn) will be required to solve the equations explicitly, which is unrealistic in the 
present computational resources. 

The net rates of lepton-number and energy exchanges between matters and 
neutrinos may not be large, and consequently an effective timescale apparently may 
not be short compared to the dynamical timescale. However, this does not immediately 
imply that one can solve the equations explicitly without bringing in any devices. For 
example the net electron capture rate vanishes in the /3-equilibrium. The achievement 
of /^-equilibrium is consequences of both cancellation of two very large weak interaction 
processes (the electron and the electron-neutrino captures) and the action of the phase 
space blocking. Note that the weak interaction processes depend enormously on the 
temperature and the lepton chemical potentials. Therefore, small error in evaluations 
of the temperature and a small deviation from the /3-equilibrium due to small error in 
estimation of the lepton fractions will produce large error and stiff source terms, and 
consequently the explicit numerical evolutions may become unstable. 

In the following of this section, I describe difficulties of implementation of weak 
interactions and neutrino cooling into the hydrodynamic equations in the conservative 
schemes in full general relativity compared with the Newtonian framework. 

In the Newtonian framework, the equations may be solved implicitly [391 HOI HU 
S21 S31 SU HH1 SHI H7] (see also [1SI E3] and references therein). The equations of 
hydrodynamics, lepton-number conservations, and neutrino processes are schematically 
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written as, 



P =0, 



Vi = S Vi (p,Y e ,T,Q„), 
Ye = Sy e (p, Y e , T, Q v ), 
e = S e (p,Y e ,T,Q v ), 
Qy = S Qv (p,Y e ,T, Q v ) 



(1) 
(2) 
(3) 
(4) 
(5) 



where p is the rest mass density, V{ is the velocity, Y e is the electron fraction, e is the 
(internal) energy of matter, T is the temperature and Q u stands for the relevant neutrino 
quantities. S"s in the right hand side stand for the relevant source terms. Comparing the 
quantities in the left-hand-side and the argument quantities in the source terms, only 
the relation between e and T is nontrivial. Usually, EOSs employed in the simulation is 
tabularized, and one dimensional search over the EOS table is required to solve them. 
To achieve this procedure in an implicit manner is quite a difficult problem. 

In the relativistic framework, the situation becomes much more complicated in 
conservative schemes, since the Lorentz factor (r) is coupled with rest mass density 
and the energy density (see Eqs. ( 1401 and ( 1471) ). and since the specific enthalpy 
(h = h(p,Y e ,T)) is coupled with the momentum (see Eq. (|45|) ). 

It should be addressed that the previous fully general relativistic works in the 
spherical symmetry [201 [21] are based on the so-called Misner-Sharp coordinates |49j. 
There are no such complicated couplings in this Lagrangian coordinates. Accordingly 
the equations may be solved essentially in the same manner as in the Newtonian 
framework. In multidimensional case, on the other hand, no Lagrangian coordinates 
are known, and the Eulerian coordinates are adopted. In the Eulerian coordinate, the 
complicated couplings inevitably appear in multidimensional case. 

Omitting the factors associated with the geometric variables (which are known 
when solving hydrodynamics equations), the equations to be solved in the relativistic 
framework are schematically written as, 



where p* is a weighted density, u a is a weighted four velocity, e is a weighted energy 
density (see Sec. 13.41 for the definition of these variables). Note that the Lorentz factor 
is appeared in the equations. 

The Lorentz factor is calculated by solving the normalization condition u a u a = — 1, 
which is rather complicated nonlinear equation schematically written as 




p*(p,T) = 0, 




(6) 
(7) 
(8) 
(9) 
(10) 



normalization [Uij t ) /normalization (^ij Pi J- ei J- 5 
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(For a solution of the equation in practice, see Sec. 13.51 ) The accurate calculation of 
the Lorentz factor and the accurate solution of the normalization condition are very 
important in the numerical relativistic hydrodynamics. 

Now, it is obvious that the argument quantities in the source terms are not simply 
related with the left-hand-side evolved quantities in Eqs. ()6l)-(TTTT). To solve the 
equations implicitly is quite difficult and there are even no successful formulations. 
Moreover it is not clear whether a convergent solution can be stably obtained numerically 
or not, since they are simultaneous nonlinear equations. Therefore, it may not be a poor 
choice to adopt an alternative approach in which the equations are solved explicitly with 
some approximations (see Sec. [3]). 

The second minor problem which is not exist in the Newtonian framework is that 
one cannot add any microphysical processes in forms of 'cooling' or 'heating' terms Q a , 
into the right hand side of hydrodynamic equations as 

V a T| = Q p . (12) 

Instead, the energy-momentum tensor of neutrinos should be introduced in the general 
relativistic framework: 

(T to %p = (T F U + (T"U, (13) 

where T tot , T F , and T v are the total energy-momentum tensor and energy-momentum 
tensor of fluid and neutrino parts, respectively (see Sec. [3] for the definition). Now, one 
should solve the following coupled equations 

V a (T F )^ = Q P , (14) 
V a (T")£ = -Qp. (15) 

Here, the source term Q a can be regarded as the local production of neutrinos through 
the weak processes, ignoring non-local neutrino capture on matter. 



3. General relativistic neutrino leakage scheme 

In the following, I describe in some detail a method for solving all of the equation in an 
explicit manner. As described in the previous section, since t wp <C tdyn in the hot dense 
matter regions, the source terms in the equations become too stiff "to be solved explicitly. 
The characteristic timescale of leakage of neutrinos from the system £i ea k> however, is 
much longer than t wp in the hot dense matter region. Note that ti ea k ~ L/c ~ t dyn where 
L is the characteristic length scale of the system. On the other hand, ti ea k is comparable 
to t wp in the free-streaming regions where the WP timescale is longer than or comparable 
with the dynamical timescale. Utilizing these facts, I approximate the original matter 
equations and reformulate them so that the source terms are to be characterized by the 
leakage timescale ti ea k- 
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3.1. Decomposition of neutrino energy-momentum tensor 

Now, the problem is that the source term Q a in Eqs. ( !T4|) and ( TTBT) becomes too stiff 'to 
solve explicitly in hot dense matter regions where t wp <C td yn - To overcome the situation, 
the following procedures are adopted. 

First, it is assumed that the energy-momentum tensor of neutrinos are be 
decomposed into 'trapped-neutrino' ((T^ T ) a ^) and 'streaming-neutrino' ((T^ s ) Q/3 ) parts 
as [27], 

{T V U = (T U ' T U + (T U ' S U- (16) 

Here, the trapped-neutrinos phenomenologically represent neutrinos which interact 
sufficiently frequently with matter and are thermalized, while the streaming-neutrino 
part describes a phenomenological flow of neutrinos streaming out of the system [27] 
(see also [38] where a more sophisticate method based on the distribution function is 
adopted in the Newtonian framework). 

Second, the locally produced neutrinos are assumed to leak out to be the streaming- 
neutrinos with a leakage rate Q„ : 

Vp{T'> a )Z = Q? k . (17) 
Then, the equation of the trapped-neutrino part becomes 

V^f a = Q a -Q^. (18) 
Third, the trapped-neutrino part is combined with the fluid part to give 

T aP = (T F ) Q/3 + CT< T ) Q/3 , (19) 
and Eqs. ([HI) and (|T8l) are combined to give 

= -QL eak - (20) 

Thus the equations to be solved is changed from Eqs. ffTll) and f[T5|) to Eqs. (T20l) 
and ( 1T71) . Note that the new equations only include the source terms Q£ which is 
characterized by the leakage timescale ti ca k- Definition of Q£ will be found in Sec. 13.31 
The energy- momentum tensor of the fluid and trapped-neutrino parts (T a p) is 
treated as that of the perfect fluid 

T a p = {p + pe + P)u a up + Pg aP , (21) 

where p and u a are the rest mass density and the 4-velocity. The specific internal energy 
density (e) and the pressure (P) are the sum of the contributions from the baryons (free 
protons, free neutrons, a-particles, and heavy nuclei), leptons (electrons, positrons, and 
trapped-neutrinos) , and the radiation as, 

P = P B + P e + P v + P r , (22) 
e = e B + e e + e v + e r , (23) 

where subscripts '£>', 'e', V, and V denote the components of the baryons, electrons 
and positrons, radiation, and trapped-neutrinos, respectively. 
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The streaming-neutrino part, on the other hand, is set to be a general form of 

(T"' s ) a /3 = En a n p + F a n p + F p n a + P aP , (24) 

where F a n a = P a pn a = 0. In order to close the system, we need an explicit expression 
of P a p. In this paper, I adopt a rather simple form P a p = x^lafi with x = 1/3. This 
approximation may work well in high density regions but will violate in low density 
regions. However, the violation will not affect the dynamics since the total amount of 
streaming-neutrinos emitted in low density regions will be small. Of course, a more 
sophisticated treatment will be necessary in a future study. 



3.2. The lepton-number conservation equations 

The conservation equations of the lepton fractions can be written schematically as 
dY R 



dt " 7e ' (25) 



dY 
dt 

dYpe 

dt 
dY vx 

dt 



>e, (26) 
7p e , (27) 
lux, (28) 



where Y e , Y ve , Y Pe , and Y vx denote the electron fraction, the electron neutrino fraction, 
the electron anti-neutrino fraction, and fi and r neutrino and anti-neutrino fractions, 
respectively. It should be addressed that, in the previous simulations based on the 
leakage schemes [2E1 EHJ I32J |37], the neutrino fractions are not solved. 

The source terms of neutrino fractions can be written, on the basis of the present 
leakage scheme, as 

Jue = I 1 ™* 1 ~ ll7\ (29) 
Ive = 7p°e Cal - TP?", (30) 

lv* = liT 1 ~ 7^ ak , (31) 

where ^ ocal and ^ eak are the local production and the leakage rates of neutrinos, 
respectively (see Sec. I3.3p . Note that only the trapped-neutrinos are responsible for 
the neutrino fractions. The thermodynamical quantities (e.g., the pressure and the 
chemical potentials) of neutrinos can be calculated from the neutrino fractions on the 
assumption of thermalization of the trapped neutrinos. 

The source term for the electron fraction conservation can be written 

7e = 7^°e cal - 7^ cal - (32) 

Since ■y l ° cal s are characterized by by the WP timescale i wp , some procedures are 
necessary to solve the lepton conservation equations explicitly. The following simple 
procedures are proposed to solve the equation stably. 



General relativistic leakage scheme 



10 



First, in each timestep n, the conservation equation of the total lepton fraction 

It = ~ lh (33) 
is solved together with the conservation equation of Y ux , Eq. (l28p. in advance of solving 
whole of the lepton conservation equations (Eqs. ( 125]) - (|28|) ). Note that the source 
term 7; = 7jf e ak — 7p e ak is characterized by the leakage timescale t\ ea ^ so that this equation 
may be solved explicitly in the hydrodynamic timescale. Then, assuming that the /3- 
equilibrium is achieved, values of the lepton fractions in the /9-equilibrium (Y& , Y@ e , and 
Yp e ) are calculated from evolved YJ. 

Second, regarding Y^ e and Yp e as the maximum allowed values of the neutrino 
fractions in the next timestep n + 1, the source terms are limited so that Y u 's in the 
timestep n+ 1 do not exceed Yf's. Then, the whole of the lepton conservation equations 
(Eqs. ff25l) - f )28|) ) are solved explicitly utilizing the limiters. 

Third, the following conditions are checked 

/J p + fie < /J n + ll ve , (34) 
f"n /^e < f^p "I" t^uei (35) 

where fi p , fi n , fi e , \i ve and fip e are the chemical potentials of protons, neutrons, electrons, 
electron neutrinos, and electron anti-neutrinos, respectively. If both conditions are 
satisfied, the values of the lepton fractions in the timestep n + 1 is set to be those in the 
/3-equilibrium value; Y@, Y@ e , and Yp e . On the other hand, if either or both conditions 
are not satisfied, the lepton fractions in the timestep n + 1 is set to be those obtained 
by solving whole of the lepton-number conservation equations. 

A limiter for the evolution of Y vx may be also necessary in some case where the pair 
processes are dominant, for example, in simulations of collapse of population III stellar 
core. In this case, the value of Y vx at the pair equilibrium (i.e. at \i vx = 0), Y^ x ir may 
be used to limit the source term. 

In the present implementation it is not necessary to somewhat artificially divide 
the system into neutrino-trapped and free-streaming regions. Therefore there is no 
discontinuous boundary which existed in the previous leakage schemes [281 1211 [32] . 

I found that simulations of the collapse of population III stellar core and the 
formation of a black hole, in which very high temperatures (T > 100 MeV) are achieved, 
can be stably performed using the simple procedure presented in this paper. 

3.3. Definition of leakage rates 

In this subsection the definitions of the leakage rates an d 7jf ak are presented. 

Because Qj, eak may be regarded as the emissivity of neutrinos measured in the fluid rest 
frame, Q 1 ^ is defined as [50] 

gteak = glcak MQ _ (36) 
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Note that although there may be a freedom to include terms H a which satisfies 
HaU a = 0, Eq. f)3.3p may be the best choice in the present framework. 

The leakage rates Ql, eak and 7^ eak are assumed to satisfy the following properties. 

(i) The leakage rates approach the local rates Q l ° cal and ^ ocal in the low density, 
transparent region. 

(ii) The leakage rates approach the diffusion rates Q^ lfr and 7^ iff in the high density, 
opaque region. 

(iii) The above two limits should be connected smoothly. 

Here, the local rates can be calculated based on the theory of weak interactions (see Sec. 
14.31 for the local rates adopted in this paper) and the diffusion rates can be determined 
based on the diffusion theory (see Sec. I4.4l for the definition of the diffusion rate adopted 
in this paper). There will be several prescriptions to satisfy the requirement (iii) [371135]. 
In this paper, the leakage rates are defined as 

gleak = (1 _ e -6r„)Qdiff + g-^Qlocd (3?) 

7 ' eak = (1 - e- bT ")^ iS + e- feTl ^ ocal , (38) 

where r u is the optical depth of neutrinos and b is a parameter which is typically set as 
b^ 1 = 2/3. The optical depth can be computed from the cross sections in a standard 
manner [371 EH] • 

3.4- Explicit forms of basic equations in leakage scheme 

The basic equations for the general relativistic hydrodynamics are the continuity 
equation, the lepton-number conservation equations, and the local conservation equation 
of the energy-momentum. The explicit forms of the equations are presented in this 
subsection for the purpose of convenience. 

3.4-1. The Continuity and lepton-number conservation equations The continuity 
equation is 

V a {pu a ) = 0. (39) 

As fundamental variables for numerical simulations, the following quantities are 
introduced: p* = pwe 6 ^ and v % = % where w = au l . Then, the continuity equation is 
written as 

d t {p*Vv) + d k {p*v k Vv) = 0, (40) 



where y/rj = ^/detr^j is the volume element of the flat space in the curvilinear 
coordinates. 

The lepton-number conservation equations (125]) - ( 12 8j) can be abbreviated as 

-w = ^ (41) 

Using the continuity equation, they become 

dt(p.Y L ^j) + d k {p*Y L v k ^q) = p* lL . (42) 



General relativistic leakage scheme 



12 



3.4-2. Energy-momentum conservation As discussed in Sec. 13.11 I solve the following 
equations. 

= -Q l ™\ (43) 



u,S\/3 



Q 



leak 

a i 



(44) 



where T a/3 and (T v,s ) a ^ are given by Eqs. fl2TT) and fl24"|) . respectively. The source term 
Qjf k is defined by Eq. (ETJ). 

As fundamental variables for numerical simulations, I define the quantities Ui = hui 
and e = hw — P(pw)~ 1 . Then, the Euler equation (7fV J gT^ ! = 7fQ Q ), and the energy 
equation (n^V/jT^ = n a Q a ) can be written as 



d t (p*u Ay /rj) + dk [p*u A v k + 



whdACt — UidA$ % + 



ae 



-I* 



■UkUidAj 



~M 



2ah(w 2 - 1) 



d, 



+ P^(«e' 



6<M 



Pae^6j 64> 



W 



dt{p*e^q) + +d k [{p*v k e + Pe e ^(v k + P k )} 

= ae^J^PK + ■^-u k u l K kl - p^ua'Wja + ae^Q a n a , 

where the subscript A denotes w or z component. 

The evolution equations of streaming- neutrinos (E and Fi) are written as 



d t { VlE) + d k [ji{aF k - (3 k E) = ^7 [aP kl K kl - F k d k a + aQ^ k n a 
d t {^F t ) + d k \^{aP k -^ k F^ 



a 



^7 ( -Ed t a + F k d^ k + 2 P dilki + ®Q 



leak 



(45) 
(46) 

(47) 

(48) 
.(49) 



3. 5. Recover of primitive variables 

In each numerical timestep, the so-called primitive variables (p, Y^, T, and Vi) and the 



Lorentz factor w = au l — Jl + ^UiUj must be calculated from the conserved quantities 
(p*, p*Yl, e, and Ui), where Yl is the representative of the lepton fractions. Since the 
equation of state (EOS) of the nuclear matter are usually tabularized in terms of the 
argument quantities (p, Y p {= Y e ), T), I am devoted to the cases of the tabularized EOS 
in the following. 

In the case where the whole of the lepton-number conservation equations are solved 
(see Sec. I3.2[) . the argument quantities (p, Y e , T) are calculated from the conserved 
quantities as follows. 

(i) Give a trial value, w, of the Lorentz factor. Then, one obtains a trial value, p, of 
the rest mass density: p = p* / {we^) . 

(ii) A trial value, T, of the temperature can be obtained by solving 

e = e E os(p, Y e , T, Y ve , Y Pe , Y ux ), (50) 
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where esos is constructed from EOS table. Note that e and esos i n general contain 
contributions from trapped-neutrinos. One dimensional search over the EOS table 
is required to obtain T. 

(iii) The next trial value of the Lorentz factor is given by solving w = 
\J\ + e~^^UiUjhr 2 , where the specific enthalpy h is calculated from EOS table 
as h = h(p, Y e , T). 

(iv) Repeat the procedures (i)-(iii) until a required degree of convergence is achieved. 
Convergent solutions of the temperature and w are obtained typically within 10 
iterations. 

On the other hand, in the case where the total lepton fraction is evolved, the 
argument quantities (p,Y e ,T) must be recovered from the conserved quantities and Y\ 
under the assumption of the /3-equilibrium. In this case, two-dimensional reconstruction 

(Yi,e) (Y e ,T) (51) 

would be required for a given w. In this case, there may be in general more than one 
couple of (Y e , T) which gives the same Y\ and e. Therefore, I adopt a different method 
to recover (p,Y e ,T) [27]. 

Under the assumption of the /3-equilibrium, the electron fraction is related to the 
total lepton fraction as Y e = Y e (p,Yi,T). Using this relation, the original EOS table 
can be reconstructed in terms of the argument quantities of (p, Yi, T). Then, the same 
strategy as in the above can be adopted. Namely, 

(i) Give a trial value, w of w. Then, one obtains a trial value, p, of the rest mass 
density. 

(ii) A trial value, T, of the temperature can be obtained by solving 

e = e E os(p, Y h f, Y vx ). (52) 
One dimensional search over the EOS table is required to obtain T. 

(iii) The next trial value of w is given by solving w = \J^ + e~ A ^^UiUjh~ 2 . 

(iv) Repeat the procedures (i)-(iii) until a required degree of convergence is achieved. 
The electron and electron neutrino fractions are given as Y e = Y e (p,Yi,T), Y ue = 
Y ue (p, Y h T), Y Pe = Y Pe (p, Y h T) in the new EOS table. 

The construction of EOS table in terms of the argument variables of (p, Yi, T) is 
important in the present implementation. 

In the case of a simplified or analytic EOS, the Newton-Raphson method may be 
applied to recover the primitive variables. In the case of tabulated EOS, by contrast, 
the Newton-Raphson method may not good approach since it requires derivatives 
of thermodynamical quantities which in general cannot be calculated precisely from 
tabulated EOS by the finite differentiating method (see also Sec. 14. 2L 
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4. Specific details of microphysics 

4-1. Equation of state 

Baryons. While any EOS table can be used in the present code, an EOS [51], [52] based 
on the relativistic mean field theory is adopted for baryon EOS (hereafter denoted by 
Shen EOS) in the present version of our code. Note that the causality is guaranteed to 
be satisfied in the relativistic EOS, whereas the sound velocity sometimes exceeds the 
velocity of the light in the non- relativistic framework, e.g., in the EOS by Lattimer and 
Swesty [53J. 

Electrons and Positrons. If a EOS table for baryons does not include the 
contributions of the leptons (electrons, positrons, and neutrinos if necessary) and 
photons, one has to consistently include these contributions to the table. Electrons 
and positrons are described as ideal Fermi gases. 

To consistently calculate the contribution of the electrons, the charge neutrality 
condition Y p = Y e should be solved in terms of the electron chemical potential p e , for 
each value of the baryon rest-mass density p and the temperature T in the EOS table: 

n e (p e , T) = n_-n + = — (53) 

m u 

in terms of p e for given p, T, and Y p . Here, m u = 931.49432 MeV is the atomic mass unit, 
and n_ and n + are the total number densities (i.e., including electron-positron pairs) of 
electrons and positrons, respectively. Then all other quantities can be calculated from 
T and p e . 

Radiations. The contribution of radiations is included in a standard manner: The 
radiation pressure and the specific internal energy density are given by 

S r = —, Pr = ^~, (54) 

P 3 

where a r is the radiation constant a r = (87r 5 fc^)(15c 3 /ip) -1 and ks and hp are the 
Boltzmann's and the Plank's constants respectively. 

Trapped- Neutrinos. In this paper, the trapped-neutrinos are assumed to interact 
sufficiently frequently with matter that be thermalized. Therefore they are described as 
ideal Fermi gases with the matter temperature. Then, from the neutrino fractions Y v , 
the chemical potentials of neutrinos are calculated by solving 

Y u = Y u (p u ,T). (55) 

Using the chemical potentials and the matter temperature, the pressure and the internal 
energy of the trapped-neutrinos are calculated. 

4-2. The sound velocity 

In high-resolution shock-capturing schemes, it is in general necessary to evaluate the 
sound velocity c s , 



c 2 = i 
s h 



dP 



dp 



P dP 

p de 



(56) 
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Here, the derivatives of the pressure are calculated by 
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dp 



dP 
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(57) 



(58) 



where , B\ V, V, and v in the sum denote contributions of the baryon, the electrons, 
radiations, and neutrinos, respectively. 

Since there are in general the phase transition regions in a EOS table for baryons and 
the EOS moreover may contain some non-smooth spiky structures, careful treatments 
are necessary when evaluating the derivatives of thermodynamical quantities. In the 
present EOS table, the derivatives are carefully evaluated so that there are no spiky 
behaviors in the resulting sound velocities. 



4-3. The local rates 

In this paper, the electron and positron captures (7^ and 7*?°) [31], the electron-positron 
pair annihilation (7^ air ) [33], the plasmon decays (7^ las ) [37], and the Bremsstrahlung 
processes (7^ rems ) [56J are considered as the local production reactions of neutrinos. 
Then, the local rates for lepton fractions are 



local 

le 


= ~ec _ pc 
lu lu 1 


(59) 


local 
' ue 


— ^ ec _|_ ^P ai r 1 ^,plas , Brems 

IU ' IU ~ IU ' IU 1 


(60) 


local 

/ ue 


- o,P c 4. o,Pair 1 ^plas , Brems 

IU ' IU 'IV 'IV ' 


(61) 


local 

/ ux 


_ ,yPair _|_ ^plas _|_ ^,Brems^ 


(62) 



Similarly, the local energy emission rate Q ° is the sum of the contributions of the 
electron and positron captures (Ql c and Q^ c ) , the electron-positron pair annihilation 
(<5£ air )) the plasmon decays (Q£ las ), and Bremsstrahlung processes (Qj^ rcms ). 



4-4- The diffusion rates 

I follow [35] for the neutrino diffusion rates 7jJ lfT and Q„ lS . I present the forms of 
the diffusion rates in the following for convenience. An alternative definition of the 
diffusion rates will be found in [3?]. The cross sections adopted in this paper are those 
of neutrino-nucleus and neutrino-nucleon scattering, and neutrino absorptions on free 
nucleons. Explicit forms of these cross sections will be found in [56J 

Ignoring the higher order corrections, the neutrino cross sections can be written in 
general as 

a{E v ) = E 2 u a, (63) 

where E u is the neutrino energy and a is 'cross section' in which E 2 dependence is 
factored out. Similarly, the opacity and the optical depth are written as 

k(E v ) = J2 ^{Eu) = El J2 Hi = E% (64) 
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T (E U ) = J n(E v )ds = E 2 V J Ms = E^f. (65) 
Now the neutrino diffusion time may be defined by [37J EH] 

2f(^) = ^^r(E u ) = Ela^t = E^ v ff\ (66) 

C Cfc 

where the distance parameter Ax(E u ) is set to be 

Ax(E u ) = a dlff ^§T. (67) 

Here a^s is a parameter which controls the diffusion rates. In this paper, I adopt 
a diff = 3 as suggested in [37j . 

Finally, the neutrino diffusion rates are defined as 

^ s Jm^=ik^ TF -M' (68) 

,diff _ f E v n v {E v ) 1 4ircg v k 2 



from which the diffusion rates Q„ and 7^ lff are easily calculated. Here, g u is the 
statistical weight factor for neutrinos and n(E v )dE v is the number density of neutrino 
in the range from E v to E v + dE v under the Fermi-Dirac distribution. 



5. Validity of general relativistic leakage scheme 

5.1. Brief summary of numerical set up 

The numerical schemes for solving the Einstein's equations are essentially same as those 
in [I]; We adopt so-called BSSN formulation [581 |59j and use 4th-order finite difference 
scheme in the spatial direction and the 3rd-order Runge-Kutta scheme in the time 
integration. The advection terms such as j3 l di(j) are evaluated by a 4th-order upwind 
scheme. The hydrodynamic equations, the lepton-number conservation equations, and 



Model 




$ c < 0.0125 


< $ c < 0.025 


< $ c < 0.05 


< $c < 0.1 


$ c >0.1 


S15 




3.26 


1.60 


0.820 


0.414 


0.217 




V 


1.005 


1.005 


1.005 


1.005 


1.005 




N 


444 


668 


924 


1212 


1532 




L (km) 


2330 


2239 


2188 


2124 


2103 


S15 


Ax 


5.10 


2.90 


1.44 


0.760 


0.396 


low 


V 


1.005 


1.005 


1.005 


1.005 


1.005 


resolution 


N 


316 


444 


636 


828 


1020 




L (km) 


2244 


2151 


2073 


2043 





Table 1. Summary of the regridding procedure. The values of the minimum grid 
spacing Axq (in units of km), the non- uniform-grid factor ry, and the grid number N 
for each range of <5> c = 1 — a c are listed. 
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Figure 1. The radial profiles of the infall velocity, the density, the entropy per baryon 
and the total lepton fraction at bounce, 2 ms and 6 ms after bounce. The results for 
the finer grid resolution (solid curve) and for the coarser grid resolution (the dotted 
curves) are shown together while they are almost identical. 

equations of streaming-neutrinos are solved using the high-resolution centered scheme 

A nonuniform grid is adopted in the numerical simulation, in which the grid spacing 
increases as 

dxj + i = rjdxj, dzi+i = rjdzi (70) 

where dxj = Xj + i — Xj, dzi = zi+i — Zi and 77 is a constant. The regridding procedure 
[60| 161] is furthermore used to compute the collapse accurately and to save the 
CPU time efficiently. For the regridding, I define an effective gravitational potential 
$ c = 1 — a c ($ c > 0) where a c is the central value of the lapse function. In Table HJ 
parameters of the regridding procedure are summarized. More detailed set up of the 
simulation will be found elsewhere [63] . 

As a test problem, I performed a collapse simulation of spherical presupernova 
core. A presupernova model (S15) of 15M with solar metallicity computed in [57] is 
adopted as the initial condition. I follow the dynamical evolution of central part which 
is composed of the Fe core and some part of the Si-shell. The density, the electron 
fraction, and the temperature are used to calculate other thermodynamical quantities 
using the EOS table. 

To check the validity of the code, the results are compared with those in the state-of- 
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X[km] 

Figure 2. Snapshots of the contours of the density (top left panels), the electron 
fraction Y e (top right panels), the entropy per baryon (bottom left panels), and the 
local neutrino energy emission rate (bottom right panels) in the x-z plane at selected 
time slices. 

the-art one-dimensional simulations (hereafter, the reference simulations) in full general 
relativity [SU IZH E5j [22] , where one dimensional general relativistic Boltzmann equation 
is solved for neutrino transfer with relevant weak interaction processes. Since neutrino 
heating processes (v e + n — > p + e~ and v e + p — > n + e + ) are not include in the present 
implementation, and multidimensional effects such as convection cannot be followed 
in the one-dimensional reference simulations, I pay particular attention in comparing 
results during the collapse and the early phase (~ 10 ms) after the bounce. After that, 
direct comparison cannot be done since in the present multidimensional code convective 
activities set in. As shown below, results in the present simulation and in the reference 
simulations agree well. 

5.2. Comparison of the radial profiles 

The collapse proceeds until the nuclear density is reached in the central part of the 
iron core. Then, the inner core experiences the bounce due to the nuclear repulsive 
forces, forming a strong shock wave at the edge of the inner core. The shock wave 
propagates outward and when it crosses the neutrino-sphere, spiky burst emissions of 
neutrinos occur (neutrino bursts): Neutrinos in hot post-shock region are copiously 
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Figure 3. Time evolution of the neutrino luminosities. The results in the finer grid 
resolution (solid curves) and in the coarser grid resolution (dashed curves) are shown 
together. The two results are almost identical until the convective phase sets in, while 
they are not in the convective phase. 



emitted without interacting matters. Eventually, negative gradients of the total lepton 
fraction are formed behind the shock since neutrinos carry away the lepton number. In 
Fig. [TJ we show the radial profiles of the infall velocity, the density, the entropy per 
baryon and the total lepton fraction at selected time slices. 

The results agree at least semi-quantitatively with those in [6U [211 EH [22] ■ In 
particular, the radial profiles of the infall velocity, the density, and the entropy per 
baryon show good agreements. No such good agreements was reported in the previous 
Newtonian simulations in which leakage schemes are adopted [211 [291 ED E21 EQ]- The 
negative gradients quantitatively are little bit steeper in the present simulation. The 
reason may be partly because the transfers of lepton-number and energy are not solved 
in the present leakage scheme. Except for this quantitative difference, the two results 
agree well. It is found that the difference can be reduced by adjusting the parameter 
aaiff introduced in Sec. 14.41 

Recall that regions of negative Y\ gradient are known to be convectively unstable 
|36j . Convective activities indeed set in in the present simulation as shown in Fig. [21 

5.3. Comparison of the Neutrino luminosities 

Comparisons of the neutrino luminosities are particularly important since they depend 
on both implementations of weak interactions (especially electron capture in the present 
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case) and treatments of neutrino cooling (the detailed leakage scheme). Also, accurate 
estimations of neutrino luminosities would be primarily important for astrophysical 
applications, since neutrinos carry away the most of energy liberated during the collapse 
as the main cooling source. 

In Fig. [31 I show neutrino luminosities calculated according to [50] 



as functions of t — bounce where bounce is time at the bounce. The result also agrees 
approximately with that in the reference simulations. The neutrino bursts occur when 
the shock wave crosses the neutrino-sphere soon after the bounce. The peak luminosity 
at the neutrino burst is L Ue « 4.5 x 10 53 ergs/s in the present simulation, which agrees 
well with that in the reference simulations. The peak luminosity and the duration 
(width) of the neutrino burst emission can be improved by adjusting the parameter 
adiff. The modulation found in the later phase t — tbounce > 10 ms is due to convective 
activities driven by negative gradients of electron fraction and entropy per baryon. 

Thus, Fig. [3] illustrates that the present detailed leakage scheme works fairly well 
and may be applied to simulations of rotating core collapse to a black hole and mergers 
of binary neutron stars. 

In the previous simulations based on the leakage scheme [321 127] where the single 
'neutrino-trapping' density is adopted, the luminosities do not agree with that in the 
reference simulations. In particular, the luminosities at the neutrino bursts are quite 
different. 

5.4- Convergence 

In Figs. CD and El I show results in the higher resolution (solid curves) and the lower 
resolution (dashed curves). The radial profiles of the two resolutions are almost identical, 
showing that convergent results are obtained in the present simulation (see Fig. [T|). In 
the time evolution of neutrino luminosities (see Fig. [3]), the two results are almost 
identical before the convective activities set in. In the later phase, on the other hand, 
the two results shows slight difference. Since the convection and the turbulence can 
occurs in a infinitesimal scale length, the smaller-scale convection and turbulence are 
captured in the finer grid resolution. Further discussions associated with the convergence 
and numerical accuracy will be found in [63]. 

6. Summary and Discussions 

6.1. Summary 

In this paper, I presented an implementation of the weak interactions and the neutrino 
cooling in the framework of full general relativity. Since the characteristic timescale of 
weak interaction processes t wp ~ \Y e /Y e \ is much shorter than the dynamical timescale 
tdyn in hot dense matters, stiff source terms appears in the equations. In general, an 




(71) 
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implicit scheme may be required to solve them [3D]- However, it is not clear whether 
implicit schemes do work or not in the relativistic framework. The Lorentz factor is 
coupled with the rest mass density and the energy density. The specific enthalpy is also 
coupled with the momentum. Due to these couplings, it is very complicated to recover 
the primitive variables and the Lorentz factor from conserved quantities. Therefore 
I proposed an explicit method to solve the equations noting that the characteristic 
timescale of neutrino leakage from the system i leak is much longer than t wp and is 
comparable to td yn - 

By decomposing the energy tensor of neutrino into the trapped-neutrino and 
the streaming-neutrino parts, the equations for the energy momentum tensor can be 
rewritten so that the source terms are characterized by the leakage timescale ti ea k (see 
Eqs. ( |20l) and (JT71)). The lepton-number conservation equations, on the other hand, 
include the source terms characterized by the WP timescale. Therefore the limiters for 
the stiff source terms are introduced to solve the lepton-number conservation equations 
explicitly (see Sec. 13 .21) . 

In the numerical relativistic hydrodynamics, it is required to calculate the primitive 
variables and the Lorentz factor from the conserved quantities. In this paper, I develop 
a robust and stable procedure for it (Sec. 13.51) . 

Finally, to check the validity of the present implementation, I performed a 
collapse simulation of spherical presupernova core and compared the results with those 
obtained in the state-of-the-art one-dimensional simulations in full general relativity 
[E31 EH [65j [22]. As shown in this paper, results in this paper agree well with those 
in the state-of-the-art simulations. Thus the present implementation will be applied to 
simulations of rotating core collapse to a black hole and mergers of binary neutron stars. 

6.2. Discussions 

Since the present implementation of the microphysics is simple and explicit, it has 
advantage that the individual microphysical processes can be easily improved and 
sophisticated. 

For example, the neutrino emission via the electron capture can be easily 
sophisticated as follows. To precisely calculate the electron capture rate, the complete 
information of the parent and daughter nuclei are required. In the nuclear equations of 
state currently available, however, a representative single-nucleus average for the true 
ensemble of heavy nuclei is adopted. The representative is usually the most abundant 
nuclei. The problem in evaluating the capture rate is that the nuclei which cause the 
largest changes in Y e are neither the most abundant nuclei nor the nuclei with the largest 
rates, but the combination of the two. In fact, the most abundant nuclei tend to have 
small rates since they are more stable than others, and the fraction of the most reactive 
nuclei tend to be small [66;, 67]. Assuming that the nuclear statistical equilibrium (NSE) 
is achieved, the electron capture rates under the NSE ensemble of heavy nuclei may be 
calculated for given (p, Y e , T). Such a numerical rate table can be easily employed in 
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the present implementation. 

Also, the neutrino cross sections can be improved. As summarized in [68] , there are 
a lot of higher order corrections to the neutrino opacities. Note that small changes in 
the opacities may result in much larger changes in the neutrino luminosities, since the 
neutrino energy emission rates strongly depend on the temperature and the temperature 
at the last scattering surface (r„ ~ crT 2 ~ 1) changes as T ~ cr -1 / 2 . Although the 
correction terms are in general very complicated, it is straightforward to include the 
corrections in the present implementation. Note that the corrections become more 
important for higher neutrino energies. Therefore, the correction terms might play roles 
in the collapse of population III stellar core and the formation of a black hole in which 
very high temperatures (T > 100 MeV) are achieved. I already started studies to explore 
the importance of these corrections in the case of black hole formation. 

As briefly described in the introduction, one of main drawbacks of the present 
implementation of the neutrino cooling is that the transfer of neutrinos are not solved. 
Although to fully solve the transfer equations of neutrinos is far beyond the scope of 
this paper, there are a lot of rooms for improvements in the treatment of the neutrino 
cooling. For example, the relativistic moment formalism [69j m particular the so- 
called Ml closure formalism, may be adopted. For this purpose, a more sophisticated 
treatment of the closure relation for P a p is required. For example, one may adopt the 
Eddington tensor of the form [7T] 



1-X , 3x-lF a F« 



: 7a/3 + 



2 w 2 F 7 Ft 



E, (72) 



where \ is the (variable) Eddington factor. In the diffusion limit where the neutrino 
pressure is isotropic x = 1/3, while in the free streaming limit x — 1- We plan to 
implement a relativistic Ml closure formalism for the neutrino transfer in the near 
future. 

To conclude, the present implementation of microphysics in full general relativity 
works fairly well. We are now in the standpoint where simulations of stellar core collapse 
to a black hole and merger of compact stellar binaries can be performed including 
microphysical processes. Fruitful scientific results will be reported in the near future. 
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